In [1]:
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import os
from shapely.geometry import Point
import geopandas as gpd
import fiona
from shapely.geometry import shape, Point
from rtree import index
import geoplot as gplt
import geoplot.crs as gcrs
import mapclassify as mc
import folium
from folium.plugins import MarkerCluster , HeatMap , HeatMapWithTime
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import seaborn as sns
import plotly.express as px
import plotly.graph_objects as go
import numpy as np
import webbrowser
In [2]:
data = pd.read_csv(os.path.join('.' , 'earthquake_data.csv') , encoding = 'utf-8')
data['地震時間'] = pd.to_datetime(data['地震時間'])
data['year'] = data['地震時間'].dt.year
data['month'] = data['地震時間'].dt.month
In [3]:
fig , ax = plt.subplots(1 , 2 , figsize = (20 , 6))
# 中小地震歷年統計
data_scale_small = data.query('2<=規模<=4.6')
data_scale_small['freq'] = None
data_scale_small = data_scale_small.groupby('year' , as_index = False).agg({'freq' : lambda x : len(x)})
sns.barplot(x = data_scale_small['year'], y = data_scale_small['freq'] , color = 'cornflowerblue' , ax = ax[0])
ax[0].set_title('Small and Medium Earthquakes')
# 嚴重地震歷年統計
data_scale_severe = data.query('規模>4.6')
data_scale_severe['freq'] = None
data_scale_severe = data_scale_severe.groupby('year' , as_index = False).agg({'freq' : lambda x : len(x)})
sns.barplot(x = data_scale_severe['year'], y = data_scale_severe['freq'] , color = 'salmon' , ax = ax[1])
ax[1].set_title('Severe Earthquakes')
Out[3]:
In [4]:
data['county'] = None
data['town'] = None
data['towncode'] = None
shape_path = os.path.join('./Taiwan_Town_Shape' , 'TOWN_MOI_1091016.shp')
town_shp = gpd.read_file(shape_path , encoding = 'utf8')
shapes = {}
properties = {}
collection = fiona.open(shape_path , encoding = 'utf-8')
for f in collection:
# town = '{},{}'.format(f['properties']['COUNTYNAME'] , f['properties']['TOWNNAME'])
towncode = int(f['properties']['TOWNCODE'])
shapes[towncode] = shape(f['geometry'])
properties[towncode] = f['properties']
# 使用R-Tree加速查詢
idx = index.Index()
for town , poly in shapes.items():
idx.insert(town , poly.bounds)
for i , (log , lat) in enumerate(data[['經度' , '緯度']].values):
print('\rProgress : {:.2%}'.format(i / len(data)) , end = '')
for towncode in idx.intersection((log , lat)):
if shapes[towncode].contains(Point(log , lat)) == True:
county = town_shp.loc[town_shp['TOWNCODE'] == '{}'.format(towncode)]['COUNTYNAME'].values[0]
town = town_shp.loc[town_shp['TOWNCODE'] == '{}'.format(towncode)]['TOWNNAME'].values[0]
data['county'].loc[i] = county
data['town'].loc[i] = town
data['towncode'].loc[i] = towncode
data.replace({None : np.nan} , inplace = True)
In [5]:
data_county = data.loc[~data['towncode'].isnull()]
data_county['freq'] = None
data_county = data_county.groupby('county' , as_index = False).agg({'freq' : lambda x : len(x)})
trace = go.Pie(labels = data_county['county'] , values = data_county['freq'] , name = 'Pie Chart' ,
hoverinfo = 'label+percent+name')
layout = go.Layout(title = 'Earthquake Freq')
fig = go.Figure(data = [trace] , layout = layout)
fig.show()
In [6]:
makercluster_map = folium.Map(location = [data['緯度'].mean(), data['經度'].mean()] , zoom_start = 7)
marker_cluster = MarkerCluster().add_to(makercluster_map)
for i in range(0 , len(data)):
if data['規模'][i] >= 4.6 :
information =\
'地震時間: {}<br>規模: {:.2f}<br>深度: {:.2f}公里<br>位置: ({} , {})'.format(data['地震時間'][i] ,
data['規模'][i] ,
data['深度'][i] ,
data['經度'][i] ,
data['緯度'][i])
m = folium.Marker(location = [data['緯度'][i] , data['經度'][i]] , popup = information)
m.add_to(marker_cluster)
makercluster_map.save(os.path.join('./web_page' , 'makercluster_map.html'))
webbrowser.open(os.path.join('./web_page' , 'makercluster_map.html'))
makercluster_map
Out[6]:
In [7]:
point_map = folium.Map(location = [data['緯度'].mean(), data['經度'].mean()] , zoom_start = 7)
for i in range(0 , len(data)):
m = folium.CircleMarker(location = [data['緯度'][i] , data['經度'][i]] ,
radius = 1 ,
color = 'red',
fill = True ,
fill_color = 'red' ,
fill_opacity = 0.4)
point_map.add_child(child = m)
point_map.save(os.path.join('./web_page' , 'point_map.html'))
webbrowser.open(os.path.join('./web_page' , 'point_map.html'))
point_map
Out[7]:
In [8]:
heat_map = folium.Map(location = [data['緯度'].mean(), data['經度'].mean()] , zoom_start = 7)
heat_map_data = []
for i in range(0 , len(data)):
heat_map_data.append([data['緯度'][i] , data['經度'][i] , float(data['規模'][i])])
heat_map.add_child(HeatMap(heat_map_data , radius = 8))
heat_map.save(os.path.join('./web_page' , 'heat_map.html'))
webbrowser.open(os.path.join('./web_page' , 'heat_map.html'))
heat_map
Out[8]:
In [9]:
data['year_month'] = data.apply(lambda x : x['地震時間'].strftime('%Y-%m') , axis = 1)
date_list = data['year_month'].sort_values().unique()
time_index = data['year_month'].sort_values().unique().tolist()
data_move = []
for year_month in date_list:
data_move.append(data.loc[data['year_month'] == year_month][['緯度' , '經度' , '規模']].values.tolist())
In [10]:
time_heatmap = folium.Map(location = [data['緯度'].mean(), data['經度'].mean()] , zoom_start = 7)
hm = HeatMapWithTime(data_move ,
index = time_index ,
radius = 10)
hm.add_to(time_heatmap)
time_heatmap.save(os.path.join('./web_page' , 'time_heatmap.html'))
webbrowser.open(os.path.join('./web_page' , 'time_heatmap.html'))
time_heatmap
Out[10]:
In [11]:
data_main_island = data.loc[~data['towncode'].isnull()]
data_main_island['towncode'] = data_main_island['towncode'].astype(int)
data_main_island['towncode'] = data_main_island['towncode'].astype(str)
poly_data = pd.merge(left = town_shp[['TOWNCODE' , 'geometry']] ,
right = data_main_island ,
left_on = 'TOWNCODE' ,
right_on = 'towncode' ,
how = 'outer')
poly_data['town with record'] = None
poly_data['town with record'].loc[poly_data['county'].isnull()] = False
poly_data['town with record'].loc[~poly_data['county'].isnull()] = True
poly_data = gpd.GeoDataFrame(poly_data)
In [12]:
point_data = poly_data.copy(deep = True)
point_data['geometry'] = [Point((0 , 0)) for _ in range(0 , len(point_data))]
for index in range(0 , len(point_data)):
point_data['geometry'].loc[index] = Point((point_data['經度'].loc[index] ,
point_data['緯度'].loc[index]))
# 簡單繪製行政區
ax = gplt.polyplot(df = town_shp ,
projection = gcrs.AlbersEqualArea() ,
edgecolor = 'white',
facecolor = 'lightgray',
linewidths = 0.5 ,
figsize = (15 , 15))
gplt.polyplot(df = poly_data.loc[point_data['town with record'] == False] ,
projection = gcrs.AlbersEqualArea() ,
hatch = '/////' , # 無地震紀錄的行政區，會在該行政區以斜線表示
edgecolor = 'white' ,
facecolor = 'gray',
ax = ax)
# 有地震紀錄的地點，所在的位置與規模會描繪於圖上
column = '規模'
cluster = 10
gplt.pointplot(df = point_data.loc[poly_data['town with record'] == True] ,
projection = gcrs.AlbersEqualArea() ,
ax = ax ,
hue = column ,
scale = column , # 以地震規模作為尺寸映射列
alpha = 0.4 , # 設置散點透明度
limits = (5 , 40) , # 設置散點的尺寸範圍
linewidths = 1 , # 散點輪廓寬度
edgecolors = 'black' ,
scheme = mc.FisherJenks((point_data.loc[point_data['town with record'] == True]['規模']).unique() , k = cluster),
cmap = 'rainbow')
# 得到mapclassify中BoxPlot的數據分層點
bp = mc.FisherJenks((point_data.loc[point_data['town with record'] == True]['規模']).unique() , k = cluster)
bins = [0] + bp.bins.tolist()
cmap = plt.get_cmap('rainbow')
# 製作圖例映射對象列表
LegendElement = [mpatches.Patch(facecolor = cmap(_ / (cluster-1)) ,
label = f'{float(bins[_])}-{float(bins[_+1])}')
for _ in range(0 , cluster)] + \
[mpatches.Patch(facecolor = 'none',
edgecolor = 'black',
linewidth = 0.2,
hatch = '*',
label = 'No Record')]
# 將制作好的圖例映射對象列表導入legend()中，並配置相關參數
ax.legend(handles = LegendElement ,
loc = 'lower right',
fontsize = 15,
title = 'Scale',
title_fontsize = 20 ,
borderpad = 0.6)
# 添加標題
plt.title('Earthquakes Scale' , fontsize = 25)
Out[12]:
